home *** CD-ROM | disk | FTP | other *** search
/ Shareware Grab Bag / Shareware Grab Bag.iso / 090 / cmln0286.arc / EXOTIC4.LTG < prev    next >
Text File  |  1986-02-03  |  2KB  |  105 lines
  1.  
  2.                             Listing 4 
  3.  
  4. define(TPI,6.2831853072)
  5. subroutine fftdif(x,y,m) #Decimation in frequency FFT--reorder output
  6. #Forward transform (negative exponential)
  7. #Arguments:
  8. #  x  Real input/output array
  9. #  y  Imaginary input/output array
  10. #  m  Dimension of x and y is 2**m
  11. dimension x(1),y(1)
  12. n=2**m
  13. kjj=n
  14. do k=1,m
  15.   {
  16.   nn=kjj
  17.   kjj=kjj/2
  18.   dthet=TPI/nn
  19.   do j=1,kjj
  20.     {
  21.     thet=(j-1)*dthet
  22.     c=cos(thet)
  23.     s=sin(thet)
  24.     do i1=j,n,nn
  25.       {
  26.       j1=i1+kjj
  27.       xs=x(i1)-x(j1)
  28.       ys=y(i1)-y(j1)
  29.       x(i1)=x(i1)+x(j1)
  30.       y(i1)=y(i1)+y(j1)
  31.       x(j1)=xs*c+ys*s
  32.       y(j1)=-xs*s+ys*c
  33.       }
  34.     }
  35.   }
  36. call reordr(x,y,m)
  37. return
  38. end
  39. subroutine fftdit(x,y,m) #Decimation in time FFT--reorder input
  40. #Back transform (Positive exponential)
  41. #Arguments:
  42. #  x  Real input/output array
  43. #  y  Imaginary input/output array
  44. #  m  Dimension of x and y is 2**m
  45. dimension x(1),y(1)
  46. call reordr(x,y,m)
  47. n=2**m
  48. nn=1
  49. do k=1,m
  50.   {
  51.   kjj=nnè  nn=nn+nn
  52.   dthet=TPI/nn
  53.   do j=1,kjj
  54.     {
  55.     thet=(j-1)*dthet
  56.     c=cos(thet)
  57.     s=sin(thet)
  58.     do i1=j,n,nn
  59.       {
  60.       j1=i1+kjj
  61.       xs=x(j1)*c-y(j1)*s
  62.       ys=x(j1)*s+y(j1)*c
  63.       x(j1)=x(i1)-xs
  64.       y(j1)=y(i1)-ys
  65.       x(i1)=x(i1)+xs
  66.       y(i1)=y(i1)+ys
  67.       }
  68.     }
  69.   }
  70. return
  71. end
  72. subroutine reordr(x,y,m) #Reorders data for FFT input or output
  73. #Arguments:
  74. #  x  Real input/output array
  75. #  y  Imaginary input/output array
  76. #  m  Dimension of x and y is 2**m
  77. dimension x(1),y(1)
  78. n=2**m
  79. do i=1,n
  80.   {
  81.   k=i-1 #Reverse-bit k to form j-1
  82.   j=1
  83.   ib=n
  84.   do l=1,m #Add bits to j-1 from top down where bits exist in k from bottom up
  85.     {
  86.     ib=ib/2
  87.     kn=k/2
  88.     j=j+(k-kn*2)*ib
  89.     k=kn
  90.     }
  91.   if(j>i) #Interchange and conjugate array elements if j>i
  92.     {
  93.     q=x(j)
  94.     x(j)=x(i)
  95.     x(i)=q
  96.     q=y(j)
  97.     y(j)=-y(i)
  98.     y(i)=-q
  99.     }
  100.   else if(i==j) #Conjugate only if j==i
  101.     y(i)=-y(i)
  102.   } #No action if j<i; elements already reordered
  103. return
  104. end
  105. è